Practical 3

st_layers(here("prac3_data", "gadm36_AUS.gpkg"))

Driver: GPKG 
Available layers:

Ausoutline <- st_read(here("prac3_data", "gadm36_AUS.gpkg"), 
                      layer='gadm36_AUS_0')

Reading layer `gadm36_AUS_0' from data source `D:\CASA\GIS\prac3_data\gadm36_AUS.gpkg' using driver `GPKG'
Simple feature collection with 1 feature and 2 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 112.9211 ymin: -55.11694 xmax: 159.1092 ymax: -9.142176
Geodetic CRS:  WGS 84

# re-project or transform CRS
AusoutlinePROJECTED <- Ausoutline %>%
  st_transform(.,3112) # GDA94, a local CRS for Australia

print(AusoutlinePROJECTED)

Simple feature collection with 1 feature and 2 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -2083066 ymin: -6460625 xmax: 2346599 ymax: -1115948
Projected CRS: GDA94 / Geoscience Australia Lambert
  GID_0    NAME_0                           geom
1   AUS Australia MULTIPOLYGON (((1775780 -64...

##From sf to sp
#AusoutlineSP <- Ausoutline %>%
#  as(., "Spatial")

##From sp to sf
#AusoutlineSF <- AusoutlineSP %>%
#  st_as_sf()

library(raster)
library(terra)
jan<-terra::rast(here("prac3_data", "wc2.1_5m_tavg_01.tif"))
# have a look at the raster layer jan
jan

class       : SpatRaster 
dimensions  : 2160, 4320, 1  (nrow, ncol, nlyr)
resolution  : 0.08333333, 0.08333333  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
source      : wc2.1_5m_tavg_01.tif 
name        : wc2.1_5m_tavg_01 
min value   :          -46.042 
max value   :           34.065 

plot(jan)

using the Mollweide projection saved to a new object. The Mollweide projection retains area proportions whilst compromising accuracy of angle and shape

# set the proj 4 to a new object

pr1 <- terra::project(jan, "+proj=moll +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs")

#or....

newproj<-"ESRI:54009"
# get the jan raster and give it the new proj4
pr1 <- jan %>%
  terra::project(., newproj)
plot(pr1)

#back to WGS 84
pr1 <- pr1 %>%
  terra::project(., "EPSG:4326")
plot(pr1)

# look in our folder, find the files that end with .tif and 
library(fs)
dir_info("D:/CASA/GIS/prac3_data/") 

# select the data we actually want
library(tidyverse)
listfiles<-dir_info("D:/CASA/GIS/prac3_data/") %>%
  filter(str_detect(path, ".tif")) %>%
  dplyr::select(path)%>%
  pull()

#have a look at the file names 
listfiles

D:/CASA/GIS/prac3_data/wc2.1_5m_tavg_01.tif D:/CASA/GIS/prac3_data/wc2.1_5m_tavg_02.tif 
D:/CASA/GIS/prac3_data/wc2.1_5m_tavg_03.tif D:/CASA/GIS/prac3_data/wc2.1_5m_tavg_04.tif 
D:/CASA/GIS/prac3_data/wc2.1_5m_tavg_05.tif D:/CASA/GIS/prac3_data/wc2.1_5m_tavg_06.tif 
D:/CASA/GIS/prac3_data/wc2.1_5m_tavg_07.tif D:/CASA/GIS/prac3_data/wc2.1_5m_tavg_08.tif 
D:/CASA/GIS/prac3_data/wc2.1_5m_tavg_09.tif D:/CASA/GIS/prac3_data/wc2.1_5m_tavg_10.tif 
D:/CASA/GIS/prac3_data/wc2.1_5m_tavg_11.tif D:/CASA/GIS/prac3_data/wc2.1_5m_tavg_12.tif 

Then load all of the data straight into a SpatRaster. A SpatRaster is a collection of raster layers with the same spatial extent and resolution.

worldclimtemp <- listfiles %>%
  terra::rast()
  
#have a look at the raster stack
worldclimtemp

class       : SpatRaster 
dimensions  : 2160, 4320, 12  (nrow, ncol, nlyr)
resolution  : 0.08333333, 0.08333333  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
sources     : wc2.1_5m_tavg_01.tif  
              wc2.1_5m_tavg_02.tif  
              wc2.1_5m_tavg_03.tif  
              ... and 9 more source(s)
names       : wc2.1~vg_01, wc2.1~vg_02, wc2.1~vg_03, wc2.1~vg_04, wc2.1~vg_05, wc2.1~vg_06, ... 
min values  :     -46.042,     -44.800,     -57.986,     -64.200,     -64.829,     -64.395, ... 
max values  :      34.065,      32.908,      33.081,      34.277,      36.299,      38.458, ... 

# access the january layer
worldclimtemp[[1]]

class       : SpatRaster 
dimensions  : 2160, 4320, 1  (nrow, ncol, nlyr)
resolution  : 0.08333333, 0.08333333  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
source      : wc2.1_5m_tavg_01.tif 
name        : wc2.1_5m_tavg_01 
min value   :          -46.042 
max value   :           34.065 

#rename the layers
month <- c("Jan", "Feb", "Mar", "Apr", "May", "Jun", 
           "Jul", "Aug", "Sep", "Oct", "Nov", "Dec")

names(worldclimtemp) <- month

# get data for January using new layer name
worldclimtemp$Jan

class       : SpatRaster 
dimensions  : 2160, 4320, 1  (nrow, ncol, nlyr)
resolution  : 0.08333333, 0.08333333  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
source      : wc2.1_5m_tavg_01.tif 
name        :     Jan 
min value   : -46.042 
max value   :  34.065 

Using a raster stack we can extract data with a single command!! For example, make a dataframe of some sample sites — Australian cities/towns.

site <- c("Brisbane", "Melbourne", "Perth", "Sydney", "Broome", "Darwin", "Orange", 
          "Bunbury", "Cairns", "Adelaide", "Gold Coast", "Canberra", "Newcastle", 
          "Wollongong", "Logan City" )
lon <- c(153.03, 144.96, 115.86, 151.21, 122.23, 130.84, 149.10, 115.64, 145.77, 
         138.6, 153.43, 149.13, 151.78, 150.89, 153.12)
lat <- c(-27.47, -37.91, -31.95, -33.87, 17.96, -12.46, -33.28, -33.33, -16.92, 
         -34.93, -28, -35.28, -32.93, -34.42, -27.64)
#Put all of this inforamtion into one list 
samples <- data.frame(site, lon, lat, row.names="site")
# Extract the data from the Rasterstack for all points 
AUcitytemp<- terra::extract(worldclimtemp, samples)

Add the city names to the rows of AUcitytemp

Aucitytemp2 <- AUcitytemp %>% 
  as_tibble()%>% 
  add_column(Site = site, .before = "Jan")

take Perth as an example. We can subset our data either using the row name:

Perthtemp <- Aucitytemp2 %>%
  filter(site=="Perth")

Make a histogram of Perth’s temperature. The tibble stored the data as double and the base hist() function needs it as numeric..

hist(as.numeric(Perthtemp))

Warning: NAs introduced by coercion

library(tidyverse)
#define where you want the breaks in the historgram
userbreak<-c(8,10,12,14,16,18,20,22,24,26)

# remove the ID and site columns
Perthtemp <- Aucitytemp2 %>%
  filter(site=="Perth")

t<-Perthtemp %>%
 dplyr::select(Jan:Dec)

  hist((as.numeric(t)), 
     breaks=userbreak, 
     col="red", 
     main="Histogram of Perth Temperature", 
     xlab="Temperature", 
     ylab="Frequency")

#Check out the histogram information R generated
histinfo <- as.numeric(t) %>%
  as.numeric()%>%
  hist(.)

histinfo

$breaks
[1] 12 14 16 18 20 22 24 26

$counts
[1] 2 2 2 1 1 2 2

$density
[1] 0.08333333 0.08333333 0.08333333 0.04166667 0.04166667 0.08333333 0.08333333

$mids
[1] 13 15 17 19 21 23 25

$xname
[1] "."

$equidist
[1] TRUE

attr(,"class")
[1] "histogram"

breaks — the cut off points for the bins (or bars), we just specified these counts — the number of cells in each bin midpoints — the middle value for each bin density — the density of data per bin

Check the layer by plotting the geometr

plot(Ausoutline$geom)

#simplify the .shp with lots of points
AusoutSIMPLE <- Ausoutline %>%
  st_simplify(., dTolerance = 1000) %>% #controls the level of generalisation in the units of the map
  st_geometry()%>%
  plot()

make sure that both of our layers are in the same coordinate reference system before combine

print(Ausoutline)

Simple feature collection with 1 feature and 2 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 112.9211 ymin: -55.11694 xmax: 159.1092 ymax: -9.142176
Geodetic CRS:  WGS 84
  GID_0    NAME_0                           geom
1   AUS Australia MULTIPOLYGON (((158.6928 -5...

#this works nicely for rasters
crs(worldclimtemp)

[1] "GEOGCRS[\"WGS 84\",\n    DATUM[\"World Geodetic System 1984\",\n        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n            LENGTHUNIT[\"metre\",1]]],\n    PRIMEM[\"Greenwich\",0,\n        ANGLEUNIT[\"degree\",0.0174532925199433]],\n    CS[ellipsoidal,2],\n        AXIS[\"geodetic latitude (Lat)\",north,\n            ORDER[1],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        AXIS[\"geodetic longitude (Lon)\",east,\n            ORDER[2],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n    ID[\"EPSG\",4326]]"

Austemp <- Ausoutline %>%
  # now crop our temp data to the extent
  terra::crop(worldclimtemp,.)

# plot the output
plot(Austemp)

exactAus<-terra::mask(Austemp, Ausoutline)
exactAus

class       : SpatRaster 
dimensions  : 551, 554, 12  (nrow, ncol, nlyr)
resolution  : 0.08333333, 0.08333333  (x, y)
extent      : 112.9167, 159.0833, -55.08333, -9.166667  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
source      : memory 
names       :     Jan,      Feb,       Mar,       Apr,    May,    Jun, ... 
min values  :  6.4875,  6.53125,  5.845714,  4.571428,  2.642, -0.394, ... 
max values  : 34.0650, 32.90800, 31.986000, 30.400000, 28.800, 27.200, ... 

#subset using the known location of the raster
hist(exactAus[[3]], col="red", main ="March temperature")

#make our raster into a data.frame to be compatible with ggplot2, using a dataframe or tibble
exactAusdf <- exactAus %>%
  as.data.frame()

library(ggplot2)
# set up the basic histogram
gghist <- ggplot(exactAusdf, 
                 aes(x=Mar)) + 
  geom_histogram(color="black", 
                 fill="white")+
  labs(title="Ggplot2 histogram of Australian March temperatures", 
       x="Temperature", 
       y="Frequency")
# add a vertical line to the hisogram showing mean tempearture
gghist + geom_vline(aes(xintercept=mean(Mar, 
                                        na.rm=TRUE)),
            color="blue", 
            linetype="dashed", 
            size=1)+
  theme(plot.title = element_text(hjust = 0.5))

put our variable (months) into a one column using pivot_longer()

squishdata<-exactAusdf%>%
  pivot_longer(
  cols = 1:12,
  names_to = "Month",
  values_to = "Temp"
)

select two month

twomonths <- squishdata %>%
  # | = OR
  filter(., Month=="Jan" | Month=="Jun")

meantwomonths <- twomonths %>%
  group_by(Month) %>%
  summarise(mean=mean(Temp, na.rm=TRUE))

meantwomonths

ggplot(twomonths, aes(x=Temp, color=Month, fill=Month)) +
  geom_histogram(position="identity", alpha=0.5)+
  geom_vline(data=meantwomonths, 
             aes(xintercept=mean, 
                 color=Month),
             linetype="dashed")+
  labs(title="Ggplot2 histogram of Australian Jan and Jun
       temperatures",
       x="Temperature",
       y="Frequency")+
  theme_classic()+
  theme(plot.title = element_text(hjust = 0.5))

data_complete_cases <- squishdata %>%
  drop_na()%>% #dropped all the NAs
  # Month column map in descending order (e.g. Jan, Feb, March..)
  mutate(Month = factor(Month, levels = c("Jan","Feb","Mar",
                                          "Apr","May","Jun",
                                          "Jul","Aug","Sep",
                                          "Oct","Nov","Dec")))

# Plot faceted histogram
ggplot(data_complete_cases, aes(x=Temp, na.rm=TRUE))+
  geom_histogram(color="black", binwidth = 5)+
  labs(title="Ggplot2 faceted histogram of Australian temperatures", 
       x="Temperature",
       y="Frequency")+
  facet_grid(Month ~ .)+
  theme(plot.title = element_text(hjust = 0.5))

# mean per month
meanofall <- squishdata %>%
  group_by(Month) %>%
  summarise(mean = mean(Temp, na.rm=TRUE))

# print the top 1
head(meanofall, n=1)

# standard deviation per month
sdofall <- squishdata %>%
  group_by(Month) %>%
  summarize(sd = sd(Temp, na.rm=TRUE))

# maximum per month
maxofall <- squishdata %>%
  group_by(Month) %>%
  summarize(max = max(Temp, na.rm=TRUE))

# minimum per month
minofall <- squishdata %>%
  group_by(Month) %>%
  summarize(min = min(Temp, na.rm=TRUE))

# Interquartlie range per month
IQRofall <- squishdata %>%
  group_by(Month) %>%
  summarize(IQR = IQR(Temp, na.rm=TRUE))

# perhaps you want to store multiple outputs in one list..
lotsofstats <- squishdata %>%
  group_by(Month) %>%
  summarize(IQR = IQR(Temp, na.rm=TRUE), 
            max=max(Temp, na.rm=T))

# or you want to know the mean (or some other stat) 
#for the whole year as opposed to each month...

meanwholeyear=squishdata %>%
  summarize(meanyear = mean(Temp, na.rm=TRUE))

Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Ctrl+Alt+I.

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Ctrl+Shift+K to preview the HTML file).

The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.